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Abstract 

We derive analytic solutions for the potential and field in a one-dimensional system of masses 
or charges with periodic boundary conditions, in other words Ewald sums for one dimension. 
We also provide a set of tools for exploring the system evolution and show that it's possible to 
construct an efficient algorithm for carrying out simulations. In the cosmological setting we show 
that two approaches for satisfying periodic boundary conditions, one overly specified and the other 
completely general, provide a nearly identical clustering evolution until the number of clusters 
becomes small, at which time the influence of any size-dependent boundary cannot be ignored. 
r^ ', Finally we compare the results with other recent work with the hope of providing clarification over 

o 

differences these issues have induced. We explain that modern formulations of physics require a 
well defined potential which is not available if the forces are screened directly. 
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I. INTRODUCTION 

One-dimensional models play an important role in physics. While they are of intrin- 
sic interest, they also provide important insights into higher-dimensional systems. One- 
dimensional plasma and gravitational systems were the first N-body systems simulated with 
early computers [1-4]. Plasma systems were used to investigate Debye screening and ther- 
modynamic equilibrium. In contrast, although the force law is similar, it was found that 
the evolution of gravitational systems toward equilibrium is extremely slow, and still is not 
completely explained [5-9]. Recently spin-offs of these models that, in some cases, are more 
amenable to computer simulation have been studied in great depth [10, 11]. Three dimen- 
sional dynamical simulations are an important component of modern cosmology. Starting 
from the precise initial conditions provided by observations of the cosmic background ra- 
diation [12], various model predictions can be compared with what we observe in "today's 
universe" [13]. It was shown by Rouet and Feix that, in common with the observed positions 
of galaxies, a cosmological version of the one-dimensional gravitational system exhibits hi- 
erarchical clustering and fractal-like behavior [14, 15]. Since then we, as well as others, have 
pursued different one- dimensional cosmological models, finding important attributes such as 
power law behavior of both the density fluctuation power spectra and two-body correlation 
function, in addition to the influence of dark energy [16-21]. In specific applications, for 
both plasma and gravity, it is preferable to adopt periodic boundary conditions because 
they avoid special treatment of the boundary and therefore best mimic a segment of the 
extended system [13, 22-24] . However, since they act as a low-pass filter, no information 
supported on wavelengths larger than the system size is available. Therefore, in the context 
of plasma and gravitational systems, as well as in the Vlasov limit [25], it is necessary to 
have a sufficiently large system that contains many Jeans' (or Debye) lengths [26]. 

In both plasma and gravitational physics, a large class of one-dimensional models are 
defined by a potential energy that satisfies Poisson's equation. In three dimensions they are 
represented by embedded systems of parallel sheets of mass or electric charge that are of 
infinite extent. From Gauss' Law, the field from such an element is directed perpendicular 
to the surface and has a constant value, independent of the distance, and proportional to the 
mass or charge density (per unit area). In a linear array of sheets of equal mass or charge 
density, the force on a given sheet is then simply proportional to the difference between 



the number of sheets on the left and right. Problems with this formulation arise in specific 
applications where it is necessary to assume periodic boundary conditions, where the motion 
takes place on a torus. A particular case arises in constructing a 1+1 dimensional model of 
the expanding universe. In the cosmological setting, astrophysicists consider a segment of 
the universe that is following the average expansion rate. They assume it is large enough to 
contain many clusters, but small enough that Newtonian dynamics is adequate to describe 
the evolution [27]. Comoving coordinates, in which the average density remains fixed, are 
employed, and it is assumed that the system obeys periodic boundary conditions [22, 23]. In 
comoving coordinates fictitious forces appear, analogous to the Coriolis force in a rotating 
system. The source of the apparent gravitational field arises from the difference between 
the actual matter distribution and a negative background density. To compute the force 
on a particle in three-dimensional systems it is possible to carry out "Ewald" sums over 
the positions of all the other particles in both the system and all its periodic "replicas" 
[24, 28]. In one dimension, because the field induced by a particle (mass or charge sheet) 
is constant, the solution is not obvious. At first glance it appears that the force is due to 
the difference between infinities and, at second glance, in the cosmological setting, that the 
negative background should exactly cancel with the force due to the particles in each replica. 
Thus the problem is fraught with ambiguity. To get a different twist, consider the fact that, 
in a periodic system, we are describing motion on a torus, so there is no clear distinction 
between left and right. In order to provide conclusive solutions to these problems, in the 
following pages we will follow the approach used by Kiessling in showing that the "so-called" 
Jeans swindle is, in fact, legitimate and not a swindle at all [29]. For clarity we will focus 
on the gravitational example. 

II. PERIODIC BOUNDARY CONDITIONS 

Let p(x) be a periodic mass density (mass density per unit length) of the one- dimensional 
system with period 1L so that g(x + 2L) = p(x). In addition write p(x) = Po + cr(x) where 
Qo represents the average of p(x) over one period and a(x) is the periodic fluctuation. If, 
instead, the total mass were bounded, we could compute the potential in the usual way. The 
potential function for a unit mass located at x' is simply 2ttG \x — x'\ so we would normally 



write 
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where the integration is over the whole line. Clearly, in the present situation, this won't 
converge. Following Kiessling [29] we introduce the screening function exp (— k \x — x'\) and 
define \l/ by 

ty(x, k) — 2-kG / \x — x'\ exp (— k \x — x'\) p(x')dx'. (2) 

For mass densities of physical interest, for example for both bounded functions as well as 
delta functions, \I/ clearly exists in the current case so long as k > . The contribution from 
the average, or background, term is quickly determined: 

d f 

^o = 2nG(——) / exp (— k \x — x'\) podx' = AnGpo/n 2 (3) 

and is independent of position. Therefore, although it blows up in the limit k — > 0, since po 
is translation invariant, ^o makes no contribution to the gravitational field. On the other 
hand we can formulate the contribution to \1/ from a(x), say \IZo-, and show that it is well 
behaved in this limit. 



A. Fourier Representation of the Potential and Field 

Since the density fluctuation o is a mass-neutral periodic function, we may represent o(x) 
as a Fourier series 

°{ X ) = E„ C n eX P(^ 7rna: /^) ( 4 ) 

where the prime indicates that there is no contribution from n — 0. Inserting into Eq. (2) 
we find 

ty ff {x, k) = 2nG y^ n c n b n exp(iirnx / L) (5) 

where 



| a; — x'\ exp (— k \x — x'\) exp(inn(x' — x)/L)dx' 
\u\ exp (—k \u\) exp(imru) / L)du 
\u\ exp (— k \u\) cos(niru/L)du 
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Now, taking the limit k — > 0, we find a general expression for the periodic potential 0(x), 

0(x) = — 47rG^ n c n (L/irn) 2 exp(innx/L) (6) 

and, for the gravitational field E(x), 

E{x) = —t— = 4irGi}2 n c n (L/irn) exp(iimx / L) . (7) 

Thus, for a periodic distribution of mass, the gravitational field is well defined. We can think 
of it as arising from both the mass in the primitive cell (— L < x < L) and the contribution 
from the infinite set of replicas or images. We see immediately that it is a solution of 
the Poisson equation, from which the result can also be obtained more directly from linear 
independence. For the sake of comparison, it is worth calculating the field contributed by 
the primitive cell alone. This is simply proportional to the difference in the mass on the 
right and left of the position x: 

E p (x) = 2ttG J a(x')[Q(x - x) - Q(x - x')} (8) 

where 6 is the usual step function. Substituting for a(x) we find 

E p (x) = AnGi^2' n c n (L/7rn)[exp(i7vnx/L) - (-1)"]. (9) 

While E p (x) and E(x) are very similar, there is an important difference: the former is forced 
to vanish at the endpoints, x = ±L for all allowed mass distributions within the primitive 
cell. This result is expected since we only take into account the field of a neutral slice. 
Therefore, in the general case, E p (x) cannot represent the field on a circle (1-torus). 

So far our treatment is quite general and applies to any one-dimensional periodic mass 
(or charge) distribution. To gain further insight and make contact with recent work, let's 
focus on the situation where the sources are 2N discrete equal-mass points (sheets) with 
positions Xj that live on the torus with the coordinate boundary points at x = L and —L 
identified. Then, in the primitive cell, the density fluctuation is 

2N 

(10) 






a p {x) = m > 
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from which we may easily calculate the Fourier coefficients, 



1 f L m 2N 

c n = — / exp (—innx'/L) a p (x')dx' = — > ^exp (—innXj/L) , (11) 

J ~L j=1 

for riy^ 0. Then the gravitational potential and field consists of the contribution from the 
primitive cell and all the replicas. From Eqs. (6,7) they reduce to 

IN oo 

cf)(x) = -AixmGL ^2 ^2 (V 7m ) 2c05 ( 7m ( a; - Xj)/L), (12) 

j=l n=l 

„, 2N oo 

E(x) = —^ = -AmG^2^2{l/n)sin{>Kn{x-x j )/L). (13) 

j=\ n=l 

Thus the Fourier representations of the periodic potential and field are straightforward. 



B. Direct Summation over Replicas 

In the case of three dimensions one cannot do much better than Eqs. (12,13) since the 
Ewald sums cannot be represented as simple analytic functions [28]. Fortunately, in the 
present case, we can improve on this situation. Consider a single particle of mass m located 
at X\. By summing over replicas, we can compute its contribution to the screened potential 
ty(x, k) directly: 

oo 

ty(x,K) = 2imiG 2_, \ x ~ x i ~~ 2rL\ exp (— k \x — x\ — 2rL\) (14) 

r=— oo 

d °° 
= 27imG(-— ) Y^ exp(-K\y 1 -2rL\) (15) 

r=—oo 

= 2nmG(-— ) V" {exp (-n(y l - 2rL)) (y± - 2rL) + exp (n(y l - 2rL)) (-y x + 2rL)} 

OK 

r=—oo 

(16) 
where y\ — x — X\ . Choose integers r<(yi)and r>(yi) such that r< < y\/2L < r> = r< + 1, 
i.e. y\/2L is bounded from below and above by this pair of adjacent integers. Then 

oo 

y exp {—n{yi — 2rL)) (yi — 2rL) = exp (— n(yi — 2r < L)) \^ exp{n2Ls) (17) 

r=— oo s=— oo 

and similarly for the second sum in Eq.(16). Therefore each of the sums in Eq.(16) can be 

evaluated in terms of a geometric series to obtain the screened potential 

d 
^ a (x, k) = 2nmG{-—) {[exp (-kY<) + exp (+«y>)] / (1 - exp (-2«L)) - 1/kL} (18) 

OK 
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where Y < — y\ — 2r < L , etc., and we have subtracted the contribution from the average or 
background density, m/2L. Evaluating the derivative and then taking the limit K — >■ 0, we 
obtain the gravitational potential 0i due to a single particle at x\. 

0i(-) = -^(^ + ^< 2 ). (19) 

It is important to recognize that <pi (x) is a periodic function of its argument and can be 
evaluated anywhere on the periodic extension of the torus, i.e. on the real line. As physicists 
we are typically interested in values of x and X\ in the primitive cell, i.e. for — L < x, X\ < L 
with the points at x — ±L identified. Then we quickly find that, for y\ > 0, < y\/2L < 1 
whereas for y\ < 0, y\/2L is sandwiched between [— 1, 0) . Either way, the potential <f>i, and 
therefore the field E\, can be represented as 

(20) 

EAx) = —-77— = 2-KirnG —(x — Xi) + Q(xi — x) — Q(x — X\) . (21) 

ox \_L 

Thus, in addition to the direct contribution from the mass located at x lt there is an additional 
quadratic term in the potential and linear term in the field. Although these contributions 
are simple, care must be taken in their interpretation. They are not, respectively, equal to 
the potential and field contributed by the component of the background located between 
the point of application x and the location of the source xi, but rather twice as large! 

A number of observations are in order. First of all, it is obvious that these functions 
reproduce exactly the Fourier series derived above for the case of a single particle. Second, 
in the limit L — > 00, they reduce to the familiar results on the line. Third, they are 
strictly functions of the displacement x — X\. This is important as all points on the torus 
are equivalent: there are no special positions or intervals. Fourth, it is not necessary to 
distinguish in which direction the distance between the points x and x\ is measured. Going 
in either direction around the torus yields the same value of 4>i- Fifth, defining 0(0) = |, 
the field vanishes at both x — X\ and at x = X\ — L mod(2L), i.e. half way around the 
torus. Finally, when xi traverses the point at L and reappears at —L, or vice-versa, there 
is no change in the field at x as we would expect from the physics. 



C. Symmetry-Based Derivation 



In the above we have employed the machinery of a screening function to obtain the desired 
result. While it has all the right properties, it is worth asking if we could have obtained it 
from a simpler route. We seek a solution of Poisson's equation for a single mass located at 

Xi , 
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where, in the primitive cell, 



The general solution is 
yielding 



a = a p (x) = m 
(f>(x) = 2irmG 



4irGa(x) 
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E(x) = - 
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(22) 

(23) 
(24) 
(25) 



where b is an arbitrary constant, as before y\ — x — x±, and we have chosen the additive 
constant in the expression for such that <p(x = X\) = 0. Symmetry requires that there 
is no preferred direction on the torus. Therefore, regardless of the location of the mass at 
Xi, the average of the field in [-L, L) must vanish. Then we immediately obtain 6 = and 
the results given in subsection II B above. We could also have arrived at this conclusion by 
noting that, for the same reason, <f>(x) can only depend on the distance between x and x\. 
Finally, the requirement that 4>(L) = <f>(—L) also demands the same conclusion. The fact 
that our limiting procedure, which is based on an exponential screening function, leads to a 
unique solution of the Poisson equation increases our confidence that the choice of a different 
screening function, e.g. a Gaussian, would not result in a different potential or field. 

For the sake of comparison, and to understand the connection with other recent work, 
it's worthwhile to consider Ei p (x), the field generated solely by a p , the charge distribution 
in the primitive cell. This would be the correct field if the net contribution from each replica 
vanished. We quickly find 



Ei p (x) 



l-nrnG 



—x + 0(xi — x) — &(x — Xi) 



(26) 



Among other problems, note that it does not satisfy the symmetry requirements discussed 
above. We will return to this formulation in the ensuing discussion. 



III. N-BODY SIMULATION 

In carrying out a simulation, we need to know the field acting on each particle. Summing 
over all the contributions, the total field at x arising from the complete system of particles 
is then simply 



E(x) = 4nmG 



N 1 

—(x-x c ) + -(N R (x)-N L (x)) 



(27) 



where, in Eq. (27), x c is the center of mass of the 2N particles in [— L, L) and Nr(x), Nl(x) 
are the number of particles to the right (left) of x counted on the segment [-L, L). Since, 
from Eq. (21), the field from a single particle vanishes at its location, Eq. (27) gives the 
correct field acting on each particle, i.e. Ej = E(x = Xj). The presence of the center of 
mass in Eq. (27) means that the instantaneous field experienced by each particle depends 
on the dipole moment of the system. The dependence on the center of mass is essential 
as it insures that when a particle passes from x = L to x = —L or vice-versa, there is no 
change in the field experienced by each particle. This was recognized as a basic problem 
in simulations of the one-dimensional, single-component plasma some time ago (see [30] 
reprinted in [1]). Perhaps the correct mathematical form of the field was not known. In 
order to avoid discontinuous jumps in the field, a polarization charge was artificially induced 
at the system boundaries, and was changed whenever a particle "switched sides". In this 
way the boundaries were seen as initially neutral reservoirs of particles. When a system 
particle enters one reservoir, another particle escapes from the other, so the boundaries are 
no longer neutral. 

For systems of interest the equations of motion of the system of particles can frequently 
be cast in the form 

^ S «,,^ + „, = £(,,) (28) 

where the value of the friction constant 7 depends on the particular model [20, 21]. In the 
cosmological setting the time has been rescaled to retain the simplicity of the equations of 
motion and increases exponentially with the co moving time coordinate [21, 31]. By carefully 
summing over the index j we find that the velocity of the center of mass obeys the simple 

equation 

dx c dv c n .. 

-^ c , ^+7^ = 0. (29) 



When a particle traverses the coordinate boundary at x = L the center of mass changes 
discontinuously. However, the center of mass velocity is a smooth function of time so the first 
order equation for v c (t) can be integrated immediately: v c (t) = v c (0) exp(-jt). In particular, 
for the special case where the center of mass is initially at rest, its velocity maintains its 
initial value. On the other hand, x c (t) will change abruptly with each boundary traversal. 

In carrying out a simulation it is necessary to obtain the crossing times of adjacent 
particles. Starting at x — —L, label the particles according to their order so that X2N > 
• • • > Xj + i > Xj > ■ ■ ■ > x\ and define Zj = Xj+i — Xj , the displacement between the 
adjacent particles. Then, from Eqs. (27, 28), we find that Zj obeys 



+ 1Wj = AirniG l—Zj-1) (30) 



__ dzj dvjj , „ f N 

Wj = ~dt' ~dT 



for j = 1, . . . , 2N — 1. To complete the ring we continue in the same sense and determine 
the rate of change of the displacement between X2N and x\, that is 2L + x\ — x 2 n = z 2 n 
, and find that it too conveniently satisfies Eq.(30). Thus, by defining a positive direction, 
or orientation, on the torus, we can keep track of the all the relative positions between 
nearest-neighbor particle pairs subject to the constraints 

27V 2JV 

5> = 2L, X>; = °- (31) 

i i 

Since, from Eq.(29) we already know the velocity of the center of mass v c , we can invert the 
set {vjj}, v c , to obtain the particle velocities Vj at any time using a matrix inversion given 
by Rybicki [32]. Conceptually, since the particles have equal mass, except for labels, they 
appear to experience elastic collisions with their neighbors. When the positions of a pair of 
particles cross, the particles exchange accelerations, but the velocities are continuous. At 
such an event the labels of the particle pair are exchanged to maintain the ordering in the 
given direction. As time progresses we see that, for this completely ordered system, when 
Zj = 0, i.e. when the j th and j + 1 st particle cross, Wj changes sign. Moreover, since the 
particle labels have been exchanged, the velocities of the two neighboring pairs, Zj-\ and 
Zj+i , also change discontinuously, i.e., Wj_i — > Wj-i + Wj and Wj + \ — > Wj +i + Wj at the 
crossing time. This is all the information required to carry out a simulation. While we 
have focused on the gravitational system of equal masses, it's straightforward to extend the 
approach to the case of unequal masses, as well as to the single and two-component plasmas. 
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Below, in Fig. 1, we present a series of snapshots from two recent simulations of a 
one-dimensional model of the expanding universe in comoving coordinates in which only 
gravitational forces apply. The model used was the one introduced originally by Rouet 
and Feix [14, 15], i.e. Eq. (28) with 7 = -4=. Each simulation employs 2 12 — 1 particles, 
and identical initial conditions were drawn from a uniform waterbag configuration in the 
/x-space. The dimensionless, scaled, unit of time is expressed in terms of the Jeans' period 
and the dimensionless length is simply the number of particles [20, 21]. It's important to 
keep in mind that the scaled time is an exponentially increasing function of the comoving 
time coordinate [21, 31]. The left sequence shows the evolution under the symmetric version 
of the system as it was originally employed (see the discussion below), while the sequence 
on the right exhibits the evolution obtained with Eq.(27). In each sequence, the left hand 
column represents a histogram of positions, while the right hand column shows the positions 
in the position-velocity plane, i.e. what statistical physicists call /z-space and astrophysicists 
call phase space. 

We observe an initial shrinkage in the phase plane due to the friction constant, followed 
by the break-up of the system into many small clusters that results from the intrinsic gravita- 
tional instability [26]. As time progresses, in each case we observe the continual, self-similar 
creation of larger clusters from smaller ones. Comparing the two snapshot sequences, we 
see that initially, and for some time, they are virtually identical. Before T = 10 there is 
no discernable difference between the two runs in the location of cluster positions in the 
phase plane. Then, at T ~ 10 , we do notice a minute difference between them. Even at 
T = 12 they are remarkably similar! Finally, at T ~ 14 , we observe a noticeable differ- 
ence between the runs: in the first sequence the few remaining clusters are slowly gathering 
towards the right hand coordinate boundary whereas, in the second, they remain evenly 
dispersed throughout the system. This occurs because there is a natural bias in the original 
implementation of the model. By forcing the system to be symmetric, the motion effectively 
takes place in a fixed external potential proportional to —x 2 (as explained below, for the 
symmetric system there is a "ghost" system of particles for x < that we don't display). 
Consequently as the system loses "energy", matter naturally gathers near the right hand 
coordinate boundary, in the neighborhood of the potential minimum. In the second simula- 
tion the potential is translation invariant on the torus so there are no favored locations and 
the clusters remain uniformly dispersed. Notice that, by the time we can see any significant 

11 
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FIG. 1: Snapshots of the density and phase space at different times for a symmetric system(left) 
and a periodic system (right). The initial positions and velocities of all particles are the same for 
both simulations. 

differences, there are only about five clusters in the system so the boundaries are starting 
to play a role in the evolution. Thus the simulation is no longer representative of the larger 
system. At earlier times, while the number of clusters was large, as we would intuitively 
anticipate, the boundaries had little effect. To illustrate the working of the algorithm, in 
Fig. 2 we plot the position of the center of mass for a small system (255 particles). We can 
see how x c shifts in time. There are periods of little change, as well as periods of significant 
change, when a cluster or group encounters a coordinate boundary. 

IV. DISCUSSION AND CONCLUSIONS 

There is a long history of work on one-dimensional plasma models, both theoretical and 
computational, dating back to the 1960's[l, 2]. Some papers explicitly discuss the case 
of periodic boundary conditions. As we pointed out earlier, in simulations of the single- 
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FIG. 2: Time evolution of the center of mass for a simulation containing 255 particles. Each time 
a particle crosses the coordinate system boundary there is a corresponding variation of unity. 

component plasma Eldridge and Feix employed a polarization field at the boundaries to 
control the discontinuities in the electric field [30]. In a review paper Kunz gives an an- 
alytical expression for the potential, but there is no derivation [33]. Periodic boundary 
conditions are required for cosmological simulations [22-24]. The first one-dimensional cos- 
mological simulations were carried out by Rouet and Feix [14, 15]. To avoid the problem 
of introducing a "polarization" field at the boundaries, they assumed that the system was 
perfectly symmetric at all times, i.e. for every particle at Xj > with velocity Vj there is an 
image or "ghost" particle at position — Xj with velocity —Vj. When a particle reaches the 
coordinate boundary, the image particle is there to meet it. Thus the 2N-particle system is 
equivalent to an N-particle system with < x < L with reflecting boundaries. Notice that 
this construction forces the center-of-mass position and velocity to vanish, simplifying the 
equations of motion (see Eq.(27)). 

In other one-dimensional simulations, Aurell et al. employed open boundaries [34, 35] . In 
their studies an initially localized fluctuation inter-penetrates a quiescent region. The field 
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they employed is essentially given by Eq. (26). Gouda and Yano [17], as well as Tatekawa 
and Maeda [16], employed the Zeldovich approximation [26]. Details concerning the type 
of boundary conditions employed were not discussed in these works but, in contrast with 
Eq. (27) above, the Zeldovich approximation, as normally derived, does not depend on 
the system center of mass. Gabrielli et al. have studied the behavior of an infinite system 
of sheets perturbed from lattice positions [19]. They also employed the screening function 
introduced by Kiessling to obtain an analytical expression for the field so, in spirit, their 
work is closely related to ours. However, there is a surprising difference in the expression 
they obtained for the gravitational field for the case of periodic boundary conditions, which 
is also given by Eq. (26). Since it lacks the explicit dependence on x c , it is not translation 
invariant on the torus and there is a discontinuity in the field when a particle passes through a 
coordinate boundary at x — ±L. Consequently it doesn't represent true motion on a torus. 
While it is tempting to contemplate that the re-introduction of particles that leave from 
x = L at x = —L , and vice-versa, is adequate to guarantee periodic boundary conditions, 
this is not the case. An additional difficulty is that the field they present is self-referential, 
i.e. since Ei p (x = x\) ^ (see Eq.(26)), the field generated by a single particle will induce 
an acceleration of itself. 

The approach taken by Gabrielli et al. and the one taken here are remarkably similar, 
so it's worth trying to sort out why they produce different results. Here, following Kiessling 
[29], we have taken the usual approach of screening the gravitational potential so in Eq.(2) 
we are simply starting with a one-dimensional version of the Yukawa potential. In contrast, 
in [19], Gabrielli et al. effectively screen the field of a particle located at x\ directly by 
exp(— K \x — Xi\). It is straightforward to verify that the potential corresponding to this 
screened field is (2irmG / k)[1 — exp(—K \x — X\\)] . Then, for an arbitrary mass distribution 
p(x), the corresponding "screened" potential is 

(2-nmG/K) / [1 - exp (-« \x - x'\)]p(x')dx' . (32) 

For an extended mass distribution, for example the periodic system considered here, it is 
apparent that this won't converge. Contemporary foundations of physics, both classical and 
quantum, are based on a Lagrangian or Hamiltonian formulation in which the potential 
plays a more fundamental role than the force. A good example is Feynman's dissertation 
where he develops the path integral [36]. There Newton's laws arise from paths for which 
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the action is an extremum. To extend the current model to, say, the quantum regime, the 
availability of a clearly defined potential, such as Eq.( 20), is essential. 

In conclusion, we have derived analytic solutions for the potential and field in a one- 
dimensional system of masses or charges with periodic boundary conditions, i.e. Ewald 
sums for one-dimension. We have seen that each particle in such a system carries with it its 
own neutralizing background, without which the potential energy cannot be defined. For a 
system of particles, we have shown that the system "polarization" or center of mass must 
be explicitly included in the force law. We have also provided a set of tools for exploring 
the system evolution and have shown that it's possible to construct an efficient algorithm 
for accomplishing this. In the cosmological setting we have shown that the difference be- 
tween the choice of completely symmetric, or just periodic, boundary conditions plays an 
insignificant role in the evolution until the number of clusters becomes small, at which time 
the influence of any boundary condition will become important. Finally, we showed that 
directly screening the force, as in [19], instead of the potential leads to a divergent potential 
function for an extended system and is therefore not suitable for the preferred formulations 
of physics based on variational priniciples. In subsequent work we will explore other settings 
where boundaries play a more prominent role. 
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